{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Latex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 使用Mindquantum框架复现量子神经网络中的贫瘠高原\n",
    "\n",
    "### 原文：Barren Plateaus in Quantum Neural Network Training Landscape, Nat. Commun(2018).\n",
    "\n",
    "李子绅\n",
    "\n",
    "邮箱：lzs03333@yeah.net\n",
    "\n",
    "华为云ID：hid_b4uryzmyfxuzzn1\n",
    "\n",
    "## 1.项目介绍\n",
    "\n",
    "### 1.1 Barren Plateaus简介\n",
    "\n",
    "&emsp;&emsp;在VQA中，很多时候会采用随机初始化线路（random parameterized quantum circuits）去优化参数。而这种方式会带来一个严重的后果，就是Barren Plateaus (BP，或称贫瘠高原)。BP其实类似于深度学习中的梯度消失，具体表现为，随问题规模增长，成本函数在参数空间中表现出指数级平滑，从而导致优化算法很难进行下去，甚至根本无法计算。产生BP的原因主要有两个，一是希尔伯特空间的维数随量子比特数量指数增长，二是由于随机初始化线路的梯度估计的复杂度。\n",
    "\n",
    "&emsp;&emsp;Andrew Arrasmith等人按如下形式严格定义了BP (arXiv:2011.12245v2)：\n",
    "\n",
    "对于$\\forall \\theta_k, \\exists b>1, \\left\\langle\\partial_kC\\left(\\vec{\\theta}\\right)\\right\\rangle_\\theta=0$，有\n",
    "\n",
    "$$Var_\\theta\\left[\\partial_kC\\left(\\vec{\\theta}\\right)\\right]\\le\\mathcal{O}\\left(\\frac{1}{b^n}\\right)\\tag{1}$$\n",
    "\n",
    "其中$n$是比特数，$C$是成本函数，$\\theta_k$表示参数，角标$\\theta$表示在参数空间中求期望。$\\left\\langle\\partial_kC\\left(\\vec{\\theta}\\right)\\right\\rangle_\\theta=0$表示成本函数的梯度在平均意义下为0。在编写代码时可以利用好这一特性计算梯度方差，因为均值为0，所以梯度的方差就是梯度的平方的期望，即\n",
    "\n",
    "$$Var_\\theta\\left[\\partial_kC\\left(\\vec{\\theta}\\right)\\right]=\\left\\langle\\left[\\partial_kC\\left(\\vec{\\theta}\\right)-\\langle\\partial_kC\\left(\\vec{\\theta}\\right)\\rangle\\right]^2\\right\\rangle_\\theta=\\left\\langle\\left[\\partial_kC\\left(\\vec{\\theta}\\right)\\right]^2\\right\\rangle_\\theta\\tag{2}$$\n",
    "\n",
    "&emsp;&emsp;目前，对于BP这个问题，有很多种解决策略被发现。比较杰出的工作是搭建量子卷积神经网络 (QCNN)，由Arthur Pesah等人提出 (arXiv:2011.02966v2 )，这种方法可以把BP中指数递减的梯度方差弱化成多项式递减。\n",
    "\n",
    "### 1.2 项目概述\n",
    "\n",
    "&emsp;&emsp;本项目主要目标是以mindquantum为框架复现这篇Nat. Commun中的BP，展示梯度方差随比特数、层数的指数衰减。mindquantum处理参数化量子线路具有很大的优势，其中求解析梯度也有相对应的接口，不仅快，精度也比自己做差分来的高，这在很大程度上降低了复现难度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 复现过程\n",
    "\n",
    "### 2.1 复现思路\n",
    "\n",
    "&emsp;&emsp;主要思路是按原文搭建量子神经网络（如下图，摘自原文），图中$R_Y(\\pi/4)$门用于制备初态，后面的含参数$R_P$门以及控制Z门合计为一个单层，这个单层在后续计算中会重复多次施加在量子线路中。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\n",
    "    box-shadow: 0 2px 4px 0 rgba(34,36,38,.12),0 2px 10px 0 rgba(34,36,38,.08);\" \n",
    "    src=\"./fig1.JPG\"\n",
    "    width= \"300px\">\n",
    "    <br>\n",
    "    <div style=\"color:orange; border-bottom: 1px solid #d9d9d9;\n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图1 用于数值模拟的随机参数化量子线路（单层）</div>\n",
    "</center>\n",
    "\n",
    "&emsp;&emsp;搭建完量子线路后，使用Mindquantum框架中get expectation with grad求解析梯度。这种方式的优势是不需要做差分，直接利用parameter shift求解析梯度，速度也比较快。求出梯度后，求其模方，然后通过多次采样取平均来使结果逼近梯度模方的期望。由2式可知，梯度模方的期望就是梯度方差，至此完成一次完整的梯度方差计算。按照原文中要求，我们只关注$\\theta_{1,1}$变化的梯度，即第一层的第一个量子门上的参数的偏导数，将其作为梯度去研究。这样做的目的是简化计算，保证代码高效运行的同时较为精确地展现出BP的图像。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 代码调用\n",
    "\n",
    "&emsp;&emsp;导入源代码和必要的包："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Import Successfully\n"
     ]
    }
   ],
   "source": [
    "import src.bp as bp\n",
    "import numpy as np\n",
    "from mindquantum import *\n",
    "import matplotlib.pyplot as plt\n",
    "print('Import Successfully')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.1 梯度方差随比特数的变化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_max = 16 # maximal number of  qubits\n",
    "p = 200 # number of layers\n",
    "hamiltonian = Hamiltonian(QubitOperator('Z0 Z1'))\n",
    "xxx = np.arange(4,n_max+1) #用于绘图时候的横坐标\n",
    "yyy = np.zeros(n_max-3)        #储存梯度方差的向量\n",
    "for n in range(4, n_max+1):\n",
    "    circuit = bp.bpansatz(n,p)\n",
    "    print(\"n = \",n)\n",
    "    yyy[n-4] = bp.get_var_partial_exp(circuit, hamiltonian)\n",
    "\n",
    "#在半指数图上绘制图像\n",
    "log_var = np.log10(yyy)\n",
    "plt.plot(xxx, log_var)\n",
    "plt.title('Var(∂E) vs # of qubits')\n",
    "plt.xlabel('number of qubits')\n",
    "plt.ylabel('Lg(Var(∂E))')\n",
    "plt.savefig('result1_var_vs_qubis.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.2 梯度方差随层数变化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian = Hamiltonian(QubitOperator('Z0 Z1'))\n",
    "for n in range(2,16,2):\n",
    "    p_max = 400 # maximal number of layers\n",
    "    xxx = np.arange(20, p_max, 20) #用于绘图时候的横坐标\n",
    "    yyy = np.zeros(int((p_max-20)/20))        #储存梯度方差的向量\n",
    "    for p in range(20, p_max, 20):\n",
    "        circuit = bp.bpansatz(n,p)\n",
    "        print(\"n=%d, p = %d\", n, p)\n",
    "        yyy[int((p-20)/20)] = bp.get_var_partial_exp(circuit, hamiltonian)\n",
    "\n",
    "#在半指数图上绘制图像\n",
    "    log_var = np.log10(yyy)\n",
    "    plt.plot(xxx, log_var)\n",
    "plt.title('Var(∂E) vs # of layers')\n",
    "plt.xlabel('number of layers')\n",
    "plt.ylabel('Lg(Var(∂E))')\n",
    "plt.savefig('result2_var_vs_layers.svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.3 结果展示\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\n",
    "    box-shadow: 0 2px 4px 0 rgba(34,36,38,.12),0 2px 10px 0 rgba(34,36,38,.08);\" \n",
    "    src=\"./fig3.JPG\"\n",
    "    width= \"700px\">\n",
    "    <br>\n",
    "    <div style=\"color:orange; border-bottom: 1px solid #d9d9d9;\n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图2 复现结果，左图是梯度方差与比特数的关系，右图是梯度方差与层数的关系，不同颜色表示不同比特数，比特数取值为2到16的偶数</div>\n",
    "</center>\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\n",
    "    box-shadow: 0 2px 4px 0 rgba(34,36,38,.12),0 2px 10px 0 rgba(34,36,38,.08);\" \n",
    "    src=\"./fig2.JPG\"\n",
    "    width= \"700px\">\n",
    "    <br>\n",
    "    <div style=\"color:orange; border-bottom: 1px solid #d9d9d9;\n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图3 原文结果，左图是梯度方差与比特数的关系，右图是梯度方差与层数的关系，不同颜色表示不同比特数，比特数取值为2到24的偶数</div>\n",
    "</center>\n",
    " \n",
    "&emsp;&emsp;左图斜率的差别应该是由固定层数的选取不同导致的；右图浮现结果看起来并没有收敛完全，可以通过增加采样次数来提高精度。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.4 潜在问题\n",
    "\n",
    "&emsp;&emsp;计算梯度方差实际上非常消耗算力，为了使样品均值收敛到真实的期望值，需要很多次采样计算，并且，为了抵消参数空间扩大带来的稀释作用还要进一步增加采样次数。这使得梯度方差的估算变得相当耗时。受制于此，再加上云平台的时间限制，本次复现也没有跑到与原文相当的规模（比特数最多到16，层数最多到380）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 项目总结\n",
    "\n",
    "&emsp;&emsp;从上面的结果可以看出，在量子神经网络中，采用随机初始化策略，会在问题规模增大时产生Barren Plateaus。这警示我们在设计算法的时候要了解问题的结构，避免采用随机的初始化策略。\n",
    "\n",
    "&emsp;&emsp;除此之外，在很多变分量子算法里，梯度方差直接与算法复杂度挂钩。因此，在做bench marking的时候，BP也是个不能忽略的点。\n",
    "\n",
    "&emsp;&emsp;在本次开源活动中，我希望能够给出一个计算梯度方差的通用方法，为mindquantum代码仓贡献一份绵薄之力。我的编程水平不高，代码写得粗糙，还有很大的优化空间。但无论如何，我还是希望能借此机会，抛砖引玉，给之后的研究者提供一点便利、进而在这个方向上做出有建设性的工作。\n",
    "\n",
    "## 未来展望\n",
    "\n",
    "&emsp;&emsp;神经网络的梯度，无论是在经典还是量子的语境下，都是一个十分重要的话题。除了梯度消失，与此相反，还有另外一个困难，即梯度爆炸。梯度消失使得成本函数在参数空间中“无处可降”；梯度爆炸使得参数更新的步长不得不减小，以适应在参数空间上快速变化的成本函数。如果参数更新的步长小于量子门的误差，将直接影响到量子神经网络在NISQ（中等噪音量子设备）上的可训练性。另外，19年有研究者提出自然梯度的概念，用成本函数在希尔伯特空间中的梯度替换传统的参数空间中的梯度，这也是个很有意思的方向。因此，我认为梯度这个方面在现在依然有很多工作需要做。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.5 64-bit",
   "name": "python37564bita1f4675f27af400381d0ebfc99217385"
  },
  "language_info": {
   "name": "python",
   "version": "3.7.5-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}